/******************************************************/
/************************ Content  ********************/
/*

** Figures **
1 - Summary statistics figures 
2 - Illustration RD and RD-DID method: First stage
3 - Illustration of mental health results
4 - Illustration of physical health and labor market results

** Tables **
5 - Complier characteristics
6 - Effects on mental health outcomes
7 - Effects on physical health outcomes
8 - Effects on labor market outcomes
9 - Effects on public health insurance costs

*/
/******************************************************/

/*********************/
/********* 1 *********/
/*********************/

clear all
cd "I:\workdata\706727\Build\Data"
use psychighfreq_AllAges_YearCollapse

keep if inrange(year,2000,2007)
gen age =floor(minagec)
keep if inrange(age,0,58)

* Create mutually exclusive categories *
replace Dbamb=0 if DPHosp==1
replace Derp=0 if DPHosp==1 | Dbamb==1 
replace Dpsychiatry=0 if Derp==1 | DPHosp==1 | Dbamb==1
replace Dpsyc=0 if Dpsychiatry==1 | Derp==1 | DPHosp==1 | Dbamb==1
replace Dtherapy=0 if Dpsyc==1 | Dpsychiatry==1 | Derp==1 | DPHosp==1 | Dbamb==1

collapse (mean) Dbamb Dpsyc Dtherapy Dpsychiatry Derp DPHosp (sum) sumDbamb=Dbamb sumDpsyc=Dpsyc sumDtherapy=Dtherapy sumDpsychiatry=Dpsychiatry sumDerp=Derp sumDPHosp=DPHosp, by(age) fast

* Create stacked outcomes *
replace Dtherapy=Dtherapy
replace Dpsyc=Dpsyc+Dtherapy
replace Dpsychiatry=Dpsychiatry+Dpsyc
replace Derp=Derp+Dpsychiatry
replace Dbamb=Dbamb+Derp
replace DPHosp=DPHosp+Dbamb

foreach var in DPHosp Dbamb Derp Dpsychiatry Dpsyc Dtherapy {
replace `var'=`var'*100
}
foreach var in DPHosp Dbamb Derp Dpsychiatry Dpsyc Dtherapy {
replace `var'=. if sum`var'<=5
}

twoway (area DPHosp Dbamb Derp Dpsychiatry Dpsyc Dtherapy age), graphregion(color(white)) ///
legend(label(1 "Inpatient Psyc. Hosp.") label(2 "Outpatient Psyc. Hosp.") label(3 "ER Psyc. Hosp.") ///
 label(4 "Psychiatrist") label(5 "Psychologist") label(6 "PCP MHC") cols(1) region(color(none)) size(small) symxsize(small) ring(0) pos(11) keygap(*0.5)) ///
 xtitle("Age") ytitle("Percent") ///
xsize(5) ysize(5) yscale(titlegap(*4))

graph export I:\Workdata\706727\Figs\MentalServicesFigure.eps , replace 
graph export I:\Workdata\706727\Figs\MentalServicesFigure.png , replace width(2000)


/* Prescription drugs */
clear all
cd "I:\workdata\706727\Build\Data"
use psychighfreq_AllAges_YearCollapse

keep if inrange(year,2000,2007)
gen age=floor(minagec)
keep if inrange(age,0,58)

gen Dotherpsykmed=(Dpsykmed==1 & Danxiety==0 & Dsleep==0 & Dadep==0)
replace Danxiety=0 if Dsleep==1 | Dadep==1 
replace Dsleep=0 if  Dadep==1 

gcollapse (mean) Dotherpsykmed Danxiety Dsleep Dadep Dantipsych Dpsychostim (sum) sumDotherpsykmed=Dotherpsykmed sumDanxiety=Danxiety sumDsleep=Dsleep sumDadep=Dadep, by(age) fast

replace Danxiety=Danxiety+Dsleep+Dadep
replace Dotherpsykmed=Dotherpsykmed+Danxiety

foreach var in  Dotherpsykmed Danxiety Dadep {
replace `var'=`var'*100
}
foreach var in  Dotherpsykmed Danxiety Dadep {
replace `var'=. if sum`var'<=5
}

twoway (area Dotherpsykmed Danxiety Dadep  age), graphregion(color(white)) ///
 xtitle("Age") ytitle("Percent") ///
 legend(label(3 "Antidepressants")  ///
		label(2 "Anxiety and tranquilizers") label(1 "Other psychiatric meds") ///
		cols(1) region(color(white)) size(small) ring(0) pos(11) symxsize(small) keygap(*0.5)) ///
		xsize(5) ysize(5) yscale(titlegap(*4))

graph export I:\Workdata\706727\Figs\MentalServicesFigure_Drugs.eps , replace 
graph export I:\Workdata\706727\Figs\MentalServicesFigure_Drugs.png , replace width(2000)


/*********************/
/********* 2 *********/
/*********************/

clear all
cd "I:\workdata\706727\Build\Data"
use psychighfreq_AllAges_final

keep if inrange(year,2000,2011) | (year==2012 & month<4)

cap drop post
gen post=((year==2008 & month>3) | inrange(year,2009,2012)) 

/* Same data as collapse */
keep if inrange(agec,10,58)

/* Regression */
g byte above=(agec<38) if !mi(agec)
g byte below=(agec>=38) if !mi(agec)

g double dist=agec-38
g double zabove=above*dist
g double zbelow=dist

g double kwt=max(0,7.5-abs(dist))

replace firstpsyc=firstpsyc*100*12

/* Collapse data and create local linear trends */
gen weight=1
collapse (sum) weight (mean) firstpsyc , by(agec post) fast

/* Fit local linear regressions */
cap drop fit37b0
qui lpoly firstpsyc agec  [aw=weight] if agec>=38 &  agec!=38 & post==0,  k(tri) nograph gen(fit37b0) degree(1) at(agec) bwidth(7.5) 
qui replace fit37b0=. if agec<38 | post==1

cap drop fit37a0
qui lpoly firstpsyc agec  [aw=weight] if agec<38 & agec>18 & agec!=38 & post==0,  k(tri) nograph gen(fit37a0) degree(1) at(agec) bwidth(7.5) 
qui replace fit37a0=. if agec>38 | agec<18 | post==1

cap drop fit18b0
qui lpoly firstpsyc agec  [aw=weight] if agec<18 &  agec!=18 & post==0,  k(tri) nograph gen(fit18b0) degree(1) at(agec) bwidth(2.5) 
qui replace fit18b0=. if agec>18 | post==1

cap drop fit37b1
qui lpoly firstpsyc agec  [aw=weight] if agec>=38 &  agec!=38 & post==1,  k(tri) nograph gen(fit37b1) degree(1) at(agec) bwidth(7.5) 
qui replace fit37b1=. if agec<38 | post==0

cap drop fit37a1
qui lpoly firstpsyc agec  [aw=weight] if agec<38 & agec>18 &  agec!=38 & post==1,  k(tri) nograph gen(fit37a1) degree(1) at(agec) bwidth(7.5) 
qui replace fit37a1=. if agec>38 | agec<18 | post==0

cap drop fit18b1
qui lpoly firstpsyc agec  [aw=weight] if agec<18 &  agec!=18 & post==1,  k(tri) nograph gen(fit18b1) degree(1) at(agec) bwidth(2.5) 
qui replace fit18b1=. if agec>18 | post==0


/* Graph */
preserve
replace agec=floor(agec*2)/2

collapse (mean) firstpsyc fit* [aw=weight], by(agec post) fast

twoway (scatter firstpsyc agec if post==1, mcolor(red) msize(small) ms(circle))  ///
	  (line fit37b1 fit37a1 fit18b1 agec  if post==1, lcolor(red*1.2 red*1.2 red*1.2) lwidth(*1.5 *1.5 *1.5)) ///
	  (scatter firstpsyc agec if post==0, mcolor(gs7) msize(small) ms(circle_hollow))  ///
	  (line fit37b0 fit37a0 fit18b0 agec  if post==0, lcolor(gs5 gs5 gs5) lwidth(*1.5 *1.5 *1.5)) ///
	  ,xline(18 38, lcolor(red red)) graphregion(color(white)) text(1.25 13.5 "RD-DID 18:" "0.7 (0.05)" , box lcolor(white) fcolor(white))  ///
	  text(1.5 42 "RD 37:""0.7 (0.03)" , box lcolor(white) fcolor(white)) ///
	  legend(order(1 "Post reform" 5 "Pre reform") rows(2) ring(0) pos(1) region(color(white))) xtitle("Age") ytitle("Percent") ///
	  xscale(titlegap(*4)) yscale(titlegap(*6)) xsize(10) ysize(5)

graph export I:\workdata\706727\Figs\Final_FiguresAndTables_RDFirstStage.png, width(2000) replace
graph export I:\workdata\706727\Figs\Final_FiguresAndTables_RDFirstStage.eps,  replace

restore

/*********************/
/********* 3 *********/
/*********************/

/** RD Psychiatric hospital **/
clear all
cd "I:\workdata\706727\Build\Data"
use psychighfreq_AllAges_final

keep if inrange(year,2009,2011) | (year==2008 & month>=4) | (year==2012 & month<4)

/* Same data as collapse */
keep if inrange(agec,20,58)

/* Regression */
g byte above=(agec<38) if !mi(agec)
g byte below=(agec>=38) if !mi(agec)

g double dist=agec-38
g double zabove=above*dist
g double zbelow=dist

g double kwt=max(0,7.5-abs(dist))

gen AnxDep=(Dbamb_AnxDep==1 | Derp_AnxDep==1 | DPHosp_AnxDep==1)
gen byte firstpsykhospAnxDep=firstpsykhosp*AnxDep

replace firstpsykhospAnxDep=firstpsykhospAnxDep*100*12

/* Post */
qui reg firstpsykhospAnxDep 1.above c.zabove c.zbelow   [pw=kwt] if  agec!=38, vce(cluster pnr)

global  dires: di %10.3f _b[1.above]
global  dibase: di %10.2f _b[_cons]
global  dise: di %10.3f _se[1.above]

/* Collapse data and create local linear trends */
gen weight=1
collapse (sum) weight (mean) firstpsykhospAnxDep , by(agec) fast

/* Fit local linear regressions */
cap drop fitb
qui lpoly firstpsykhospAnxDep agec  [aw=weight] if agec>=38 &  agec!=38,  k(tri) nograph gen(fitb) degree(1) at(agec) bwidth(7.5) 
qui replace fitb=. if agec<38

cap drop fita
qui lpoly firstpsykhospAnxDep agec  [aw=weight] if agec<38 &  agec!=38,  k(tri) nograph gen(fita) degree(1) at(agec) bwidth(7.5) 
qui replace fita=. if agec>38


/* Graph */
preserve
replace agec=floor(agec*2)/2

collapse (mean) firstpsykhospAnxDep fita fitb [aw=weight], by(agec) fast

twoway (scatter firstpsykhospAnxDep agec , mcolor(red) msize(small))  ///
	  (line fita fitb agec , lcolor(red*1.2 red*1.2) lwidth(*1.5 *1.5)) ///
	  ,xline(38, lcolor(black)) graphregion(color(white)) text(0.5 28 "RD 37:""-0.04 (0.016)" , box lcolor(white) fcolor(white)) ///
	  xtitle("Age") ytitle("Percent") legend(order(1 "Post reform") holes(1 2) colfirst cols(3) ring(0) pos(1) region(color(white))) ///
	  xscale(titlegap(*4)) yscale(titlegap(*6)) xsize(5) ysize(5)

graph export I:\workdata\706727\Figs\Final_FiguresAndTables_RD_PsycHosp.png, width(2000) replace
graph export I:\workdata\706727\Figs\Final_FiguresAndTables_RD_PsycHosp.eps,  replace

restore

/** RD-DID Suicide attempts **/
clear
cd "I:\workdata\706727\Build\Data"
use psychighfreq_AllAges_final

keep if inrange(year,2000,2012)
drop if year==2012 & month>3
keep if inrange(agec,10,30) 

/* Convert to yearly rates */
foreach var in  firstsui    {
replace `var'=`var'*100*12
}

cap drop Dyear
gen Dyear=0 if inrange(year,2000,2007) | (year==2008 & month<4)
replace Dyear=1 if inrange(year,2009,2012) | (year==2008 & month>=4)

gcollapse (mean) firstsui , by(agec Dyear) fast

rename agec age
rename firstsui koef

preserve
keep if inrange(age,10,26)

gen fita0=.
gen fita1=.
gen fitb0=.
gen fitb1=.
foreach var in sui  {
	 
	foreach y in 0 1 {
		
	/* Fit local linear regressions */
	cap drop temp
	qui lpoly koef age   if age>18 &  age!=18 & Dyear==`y' ,  k(tri) nograph gen(temp) degree(1) at(age) bwidth(2.5) 
	qui replace fita`y'=temp if age>=18 & Dyear==`y' 

	cap drop temp
	qui lpoly koef age   if age<18 &  age!=18 & Dyear==`y' ,  k(tri) nograph gen(temp) degree(1) at(age) bwidth(2.5) 
	qui replace fitb`y'=temp if age<=18 & Dyear==`y'
	}
} 
replace age=floor(age*2)/2
gcollapse (mean) koef fita0 fita1 fitb0 fitb1, by(age Dyear)	fast

foreach var in sui {

twoway (scatter koef age if Dyear==0, mcolor(gs7) msymbol(circle_hollow)) (scatter koef age if Dyear==1 , mcolor(red))  ///
		(line fitb0 age if Dyear==0  & age<=18, lcolor(gs7*1.4) lwidth(*1.5)) ///
		(line fita0 age if Dyear==0  & age>=18, lcolor(gs7*1.4) lwidth(*1.5))  ///
		(line fitb1 age if Dyear==1  & age<=18, lcolor(red*1.4) lwidth(*1.5)) ///
		(line fita1 age if Dyear==1  & age>=18, lcolor(red*1.4) lwidth(*1.5))  ///		
,graphregion(color(white)) xtitle("Age") ytitle("Percent") legend(order(2 "Post reform" 1 "Pre reform") rows(2) ring(0) region(color(white) margin(zero)) size(*0.9) pos(5)) xline(18, lcolor(black)) xlabel(10(2)26)  xsize(5) ysize(5) text(0.235 14 "RD-DID 18:" "-0.06 (0.023)" , box lcolor(white) fcolor(white))

graph export  "I:\workdata\706727\Figs\Final_FiguresAndTables_RDDID_Sui.png", width(2000) replace
graph export  "I:\workdata\706727\Figs\Final_FiguresAndTables_RDDID_Sui.eps", replace

}
restore


/*********************/
/********* 4 *********/
/*********************/

/***** ES-RD *******/
clear 
cd "I:\workdata\706727\Build\Data"
use PsykTestLongRunRD, clear

keep if inrange(treatagec,18,50)
keep if inrange(distyrm,-5,7)

keep pnr dupli distyrm treatagec tyrm yrm agec inpat Dpsyc Ddisab Dsick

replace inpat=0 if inpat<0 | mi(inpat)

/* Winzorize */
qui sum inpat if inpat>0 & !mi(inpat), d
replace inpat=r(p99) if inpat>r(p99) & !mi(inpat)

egen pnrg=group(pnr dupli)

replace yrm=floor(yrm*4)/4
gen tyrm2=floor(tyrm*4)/4
cap drop distyrm*
gen distyrm=yrm-tyrm2

gcollapse (mean) inpat Dpsyc Ddisab Dsick, by(pnr pnrg distyrm treatagec tyrm yrm) fast

keep if inrange(tyrm,2008.25,2012.25)
keep if inrange(treatagec,30,47)

g byte postshock=(distyrm>=0) & !mi(distyrm)

egen distg=group(distyrm)
egen yrmg=group(yrm)

sum distg 
global max=r(max)
global maxm1=r(max)-1

/* set base group to three periods before treatment */
sum distg if distyrm==0
global bdist=3
global base=r(max)-${bdist}

cap drop bgdistg
gen bgdistg=distg
replace bgdistg=${max} if distg==${base}
replace bgdistg=${base} if distg==${max}
replace bgdistg=${max}-1 if distg==1
replace bgdistg=1 if distg==${max}-1

cap drop d above zabove zbelow 
g double d=treatagec-38
g byte above=(treatagec<38) if !mi(treatagec)
g double zabove=above*d
g double zbelow=d

cap drop kwt
g double kwt=max(0,7.5-abs(d))

** Load slightly different version of areg (only appropriate when clusters nest absorb variable) **
cap program drop aregben Estimate Replay
do "I:\Workdata\706727\STATA Programs\aregben.do"

cap postclose john
postfile john str30 var dist treat koef lower upper cons using "I:\workdata\706727\Figs\PsykTestLongRunRD_EventGraphs_RD", replace

local dmax=${max}

foreach var in  Dpsyc Dsick Ddisab inpat {

** Slightly different base group for inpatient contacts (five quarters before shock) **
if ("`var'"=="inpat") {
sum distg if distyrm==0
global bdist=5
global base=r(max)-${bdist}

cap drop bgdistg
gen bgdistg=distg
replace bgdistg=${max} if distg==${base}
replace bgdistg=${base} if distg==${max}
replace bgdistg=${max}-1 if distg==1
replace bgdistg=1 if distg==${max}-1
}


qui reg `var' 1.above c.zabove c.zbelow [pw=kwt] if treatage!=38 & distg==${base}
global cons=_b[_cons]+_b[1.above]

qui aregben `var' ibn.yrmg  ib${maxm1}.bgdistg ib${maxm1}.bgdistg#(c.zabove c.zbelow) ib${maxm1}.bgdistg#1.above [pw=kwt] if treatage!=38 ,vce(cluster pnr) absorb(pnrg)
 
mat V=e(V)
mat b=e(b)

	forvalues d=1/`dmax' {

	global cont=_b[`d'.bgdistg]
	global secont=_se[`d'.bgdistg]
	global treat=_b[`d'.bgdistg#1.above]+_b[`d'.bgdistg]
	global setreat=sqrt(V["`d'.bgdistg#1.above","`d'.bgdistg#1.above"]+V["`d'.bgdistg","`d'.bgdistg"]+2*V["`d'.bgdistg#1.above","`d'.bgdistg"])
	
	post john ("`var'") (`d') (0) (${cont}) (${cont}-1.96*${secont}) (${cont}+1.96*${secont}) (${cons})
	post john ("`var'") (`d') (1) (${treat}) (${treat}-1.96*${setreat}) (${treat}+1.96*${setreat}) (${cons})
	}
}
postclose john


/** Create graphs **/

** First stage **
clear all
cd "I:\Workdata\706727\Figs"
use PsykTestLongRunRD_EventGraphs_RD

set scheme s2color

gen olddist=dist
replace dist=18 if olddist==49
replace dist=49 if olddist==18 
replace dist=1 if olddist==48  
replace dist=48 if olddist==1  

sort var dist 

replace dist=(dist-21)/4

replace koef=koef*100 if var=="Dpsyc"
replace lower=lower*100 if var=="Dpsyc"
replace upper=upper*100 if var=="Dpsyc"
replace cons=cons*100 if var=="Dpsyc"

** Psyc **
qui sum cons if var=="Dpsyc"
global  dibase: di %10.1f r(mean)

twoway 	(rcap lower upper dist if var=="Dpsyc"  & treat==0, color(black))  ///
		(rcap lower upper dist if var=="Dpsyc"  & treat==1, color(red))  ///
		(connect koef dist if var=="Dpsyc"  & treat==0, color(black) mfcolor(white))  ///
		(connect koef dist if var=="Dpsyc"  & treat==1, color(red) ) ///
		,graphregion(color(white)) xtitle("Years from mental health shock",size(medlarge)) ytitle("Percentage points",size(medlarge)) ///
legend(order(4 "Treatment" 3 "Control") rows(2) ring(0) pos(11) region(color(white) margin(small)) size(medium) rowgap(*0.5)) ///
text(7.5 -2.75 "Baseline:${dibase}" , box lcolor(white) fcolor(white) size(medlarge)) ///
xlabel(-5(1)7,labsize(medlarge)) ylabel(,labsize(medlarge)) xscale(titlegap(*4)) xsize(5) ysize(4) xline(0, lcolor(black)) name(Dpsyc, replace)

graph export "I:\Workdata\706727\Figs\PsykTestLongRunRD_FinalIlluEventGraph_Dpsyc.png", replace width(2000)
graph export "I:\Workdata\706727\Figs\PsykTestLongRunRD_FinalIlluEventGraph_Dpsyc.eps", replace 

** inpatient hosp **
clear all
cd "I:\Workdata\706727\Figs"
use PsykTestLongRunRD_EventGraphs_RD

set scheme s2color

gen olddist=dist
replace dist=16 if olddist==49
replace dist=49 if olddist==16 
replace dist=1 if olddist==48  
replace dist=48 if olddist==1  

sort var dist 

replace dist=(dist-21)/4

sum cons if var=="inpat"
global  dibase: di %10.3f r(mean)

twoway 	(rcap lower upper dist if var=="inpat"  & treat==0, color(black))  ///
		(rcap lower upper dist if var=="inpat"  & treat==1, color(red))  ///
		(connect koef dist if var=="inpat"  & treat==0, color(black) mfcolor(white))  ///
		(connect koef dist if var=="inpat"  & treat==1, color(red) ) ///
		,graphregion(color(white)) xtitle("Years from mental health shock",size(medlarge)) ytitle("Number of services",size(medlarge)) ///
legend(order(4 "Treatment" 3 "Control") rows(2) ring(0) pos(11) region(color(white) margin(small)) size(medlarge) rowgap(*0.5)) ///
text(0.0075 -2.75 "Baseline:${dibase}" , box lcolor(white) fcolor(white) size(medlarge)) ///
xlabel(-5(1)7,labsize(medlarge)) ylabel(-0.01(0.005)0.02,labsize(medlarge)) xscale(titlegap(*4)) xsize(5) ysize(4) xline(0, lcolor(black)) name(inpat, replace)

graph export "I:\Workdata\706727\Figs\PsykTestLongRunRD_FinalIlluEventGraph_inpat.png", replace width(2000)
graph export "I:\Workdata\706727\Figs\PsykTestLongRunRD_FinalIlluEventGraph_inpat.eps", replace 


/*** Labor Market ***/
clear all
cd "I:\Workdata\706727\Figs"
use PsykTestLongRunRD_EventGraphs_RD

set scheme s2color

gen olddist=dist
replace dist=18 if olddist==49
replace dist=49 if olddist==18 
replace dist=1 if olddist==48  
replace dist=48 if olddist==1  

sort var dist 

replace dist=(dist-21)/4

/* Convert to percentage points */
foreach var in koef lower upper cons {
replace `var'=`var'*100
}

** Sickness benefits **
sum cons if var=="Dsick"
global  dibase: di %10.1f r(mean)

twoway 	(rcap lower upper dist if var=="Dsick"  & treat==0, color(black))  ///
		(rcap lower upper dist if var=="Dsick"  & treat==1, color(red))  ///
		(connect koef dist if var=="Dsick"  & treat==0, color(black) mfcolor(white))  ///
		(connect koef dist if var=="Dsick"  & treat==1, color(red) ) ///
		,graphregion(color(white)) xtitle("Years from mental health shock",size(medlarge)) ytitle("Percentage points",size(medlarge)) ///
legend(order(4 "Treatment" 3 "Control") rows(2) ring(0) pos(11) region(color(white) margin(small)) size(medlarge) rowgap(*0.5)) ///
text(5 -2.75 "Baseline:${dibase}" , box lcolor(white) fcolor(white) size(medlarge)) ///
xlabel(-5(1)7,labsize(medlarge)) ylabel(,labsize(medlarge))  xscale(titlegap(*4)) xsize(5) ysize(4) xline(0, lcolor(black)) name(inpat, replace)

graph export "I:\Workdata\706727\Figs\PsykTestLongRunRD_FinalIlluEventGraph_Dsick.png", replace width(2000)
graph export "I:\Workdata\706727\Figs\PsykTestLongRunRD_FinalIlluEventGraph_Dsick.eps", replace 

** Disability pension **
sum cons if var=="Ddisab"
global  dibase: di %10.1f r(mean)

twoway 	(rcap lower upper dist if var=="Ddisab"  & treat==0, color(black))  ///
		(rcap lower upper dist if var=="Ddisab"  & treat==1, color(red))  ///
		(connect koef dist if var=="Ddisab"  & treat==0, color(black) mfcolor(white))  ///
		(connect koef dist if var=="Ddisab"  & treat==1, color(red) ) ///
		,graphregion(color(white)) xtitle("Years from mental health shock",size(medlarge)) ytitle("Percentage points",size(medlarge)) ///
legend(order(4 "Treatment" 3 "Control") rows(2) ring(0) pos(11) region(color(white) margin(small)) size(medlarge) rowgap(*0.5)) ///
text(5 -2.75 "Baseline:${dibase}" , box lcolor(white) fcolor(white) size(medlarge)) ///
xlabel(-5(1)7,labsize(medlarge)) ylabel(,labsize(medlarge))  xscale(titlegap(*4)) xsize(5) ysize(4) xline(0, lcolor(black)) name(inpat, replace)

graph export "I:\Workdata\706727\Figs\PsykTestLongRunRD_FinalIlluEventGraph_Ddisab.png", replace width(2000)
graph export "I:\Workdata\706727\Figs\PsykTestLongRunRD_FinalIlluEventGraph_Ddisab.eps", replace 


/*********************/
/********* 5 *********/
/*********************/

/* RD - DID */
clear
cd "I:\workdata\706727\Build\Data"
use psychighfreq_AllAges_finalfull

keep if inrange(year,1998,2011)
drop if year==2011 & month>3
keep if inrange(agec,13,23) 

/* Generate variables */

gen cohab=(!mi(efalle))
gen byte famany_fam=((famany_child>0 & !mi(famany_child)) | (famany_par>0 & !mi(famany_par)) | (famany_sib>0 & !mi(famany_sib)))

gen anymental=(Dpsyc==1 | Dpsychiatry==1 | Dpsykhosp==1 | Dtherapy==1 | Dpsyk_test==1)

/* Merge with own income, education and labor market transfers one year prior */
rename year year2
gen year=year-1
merge m:1 pnr year using incdata, keepusing(pnr year inc pstill isced)
drop if _merge==2
drop _merge

merge m:1 pnr year using J:\Workdata\706727\Temp\Laborcollapse, keepusing(pnr year Dallnonparleave)
drop if _merge==2
drop year _merge
rename year2 year

gen byte EmpEduc=(pstill=="11" | pstill=="12" | pstill=="13" | pstill=="14" | pstill=="19" | pstill=="20" | pstill=="31" | pstill=="32" | pstill=="33" | pstill=="34" | pstill=="35" | pstill=="36"  | pstill=="37" | pstill=="91" )
gen byte highschool=(isced>=3 & !mi(isced))

* Use values from one year prior *
sort pnr year month
by pnr: gen lagfamany=famany[_n-12] if pnr==pnr[_n-12] 
replace famany=lagfamany

by pnr: gen lagfamany_fam=famany_fam[_n-12] if pnr==pnr[_n-12] 
replace famany_fam=lagfamany_fam

by pnr: gen lagcohab=cohab[_n-12] if pnr==pnr[_n-12] 
replace cohab=lagcohab

keep if  !mi(agec) & !mi(koen) & !mi(EmpEduc) & !mi(famany) & !mi(famany_fam) & !mi(inc)  & !mi(Dallnonparleave) & !mi(highschool) & !mi(cohab) & !mi(year)

* Merge with severity prediction *
merge m:1 pnr year using "I:/Workdata/706727/Build/Data/MLSevdata_Pred.dta"
keep if _merge==3
drop _merge

cap drop rankage
gen rankage=floor(agec)

sort year rankage pred
by year rankage: gen double predrank=_n
by year rankage: gen double maxpredrank=_N
replace predrank=ceil((((predrank-0.5)/maxpredrank)*100)/(100/4))
rename predrank highrisk

/* Create income ranks */
set seed 112
gen double random=runiform()
replace random=(random-0.5)/_N

cap drop nonmiss
gen nonmiss=(!mi(inc) & !mi(rankage) & !mi(koen))

sort nonmiss year rankage koen inc random
by nonmiss year rankage koen: gen double incrank=_n
by nonmiss year rankage koen: gen double maxincrank=_N
replace incrank=ceil((((incrank-0.5)/maxincrank)*100+random)/(100/2))
replace incrank=. if nonmiss==0

* Same sample as year_collapse *
keep if inrange(agec,15,23) & inrange(year,2000,2011)

/* Regression */
cap drop Dyear
gen Dyear=0 if inrange(year,2000,2007) | (year==2008 & month<4)
replace Dyear=1 if inrange(year,2009,2012) | (year==2008 & month>=4)

g byte above=(agec>=18) if !mi(agec)
g byte below=(agec<18) if !mi(agec)

g double dist=agec-18
g double zabove=above*dist
g double zbelow=dist

g double kwt=max(0,2.5-abs(dist))

/** overall **/
qui reg firstpsyc 1.above c.zabove c.zbelow 1.Dyear (1.above c.zabove c.zbelow)#1.Dyear   [pw=kwt] if  agec!=18
global overall=_b[1.above#1.Dyear]

qui reg anymental  1.above c.zabove c.zbelow 1.Dyear (1.above c.zabove c.zbelow)#1.Dyear   [pw=kwt] if  agec!=18
global overall_AT=_b[_cons]+_b[1.above]+_b[1.Dyear]

/** Gender **/
qui reg firstpsyc  1.above c.zabove c.zbelow 1.Dyear (1.above c.zabove c.zbelow)#1.Dyear    [pw=kwt] if  agec!=18 & koen==2
global pc_RDDID_gender_2=(_b[1.above#1.Dyear]/${overall})*100

qui reg anymental  1.above c.zabove c.zbelow 1.Dyear (1.above c.zabove c.zbelow)#1.Dyear [pw=kwt] if  agec!=18 & koen==2
global pt_RDDID_gender_2=((_b[_cons]+_b[1.above]+_b[1.Dyear])/${overall_AT})*100

/** Single **/
qui reg firstpsyc  1.above c.zabove c.zbelow 1.Dyear (1.above c.zabove c.zbelow)#1.Dyear  [pw=kwt] if  agec!=18 & cohab==0
global pc_RDDID_single=(_b[1.above#1.Dyear]/${overall})*100

qui reg anymental  1.above c.zabove c.zbelow 1.Dyear (1.above c.zabove c.zbelow)#1.Dyear [pw=kwt] if  agec!=18 & cohab==0
global pt_RDDID_single=((_b[_cons]+_b[1.above]+_b[1.Dyear])/${overall_AT})*100

/** Income Tertile **/
foreach t in 1  {

qui reg firstpsyc  1.above c.zabove c.zbelow 1.Dyear (1.above c.zabove c.zbelow)#1.Dyear  [pw=kwt] if  agec!=18 & incrank==`t'
global pt_RDDID_inc_`t'=(_b[1.above#1.Dyear]/${overall})*100

qui reg anymental  1.above c.zabove c.zbelow 1.Dyear (1.above c.zabove c.zbelow)#1.Dyear [pw=kwt] if  agec!=18 &  incrank==`t'
global pc_RDDID_inc_`t'=((_b[_cons]+_b[1.above]+_b[1.Dyear])/${overall_AT})*100
}

/** High school education **/
qui reg firstpsyc  1.above c.zabove c.zbelow 1.Dyear (1.above c.zabove c.zbelow)#1.Dyear  [pw=kwt] if  agec!=18 & highschool==1
global pc_RDDID_educ=(_b[1.above#1.Dyear]/${overall})*100

qui reg anymental 1.above c.zabove c.zbelow 1.Dyear (1.above c.zabove c.zbelow)#1.Dyear [pw=kwt] if  agec!=18 & highschool==1
global pt_RDDID_educ=((_b[_cons]+_b[1.above]+_b[1.Dyear])/${overall_AT})*100

/** Labor market transfers **/
qui reg firstpsyc 1.above c.zabove c.zbelow 1.Dyear (1.above c.zabove c.zbelow)#1.Dyear  [pw=kwt] if  agec!=18 & Dallnonparleave==1
global pc_RDDID_labortrans=(_b[1.above#1.Dyear]/${overall})*100

qui reg anymental  1.above c.zabove c.zbelow 1.Dyear (1.above c.zabove c.zbelow)#1.Dyear [pw=kwt] if  agec!=18 & Dallnonparleave==1
global pt_RDDID_labortrans=((_b[_cons]+_b[1.above]+_b[1.Dyear])/${overall_AT})*100

/** Own mental health problems **/
qui reg firstpsyc 1.above c.zabove c.zbelow 1.Dyear (1.above c.zabove c.zbelow)#1.Dyear  [pw=kwt] if  agec!=18 & famany==1
global pc_RDDID_famany=(_b[1.above#1.Dyear]/${overall})*100

qui reg anymental  1.above c.zabove c.zbelow 1.Dyear (1.above c.zabove c.zbelow)#1.Dyear [pw=kwt] if  agec!=18 & famany==1
global pt_RDDID_famany=((_b[_cons]+_b[1.above]+_b[1.Dyear])/${overall_AT})*100

/** Own mental health problems **/
qui reg firstpsyc 1.above c.zabove c.zbelow 1.Dyear (1.above c.zabove c.zbelow)#1.Dyear  [pw=kwt] if  agec!=18 & famany_fam==1
global pc_RDDID_famany_fam=(_b[1.above#1.Dyear]/${overall})*100

qui reg anymental  1.above c.zabove c.zbelow 1.Dyear (1.above c.zabove c.zbelow)#1.Dyear [pw=kwt] if  agec!=18 & famany_fam==1
global pt_RDDID_famany_fam=((_b[_cons]+_b[1.above]+_b[1.Dyear])/${overall_AT})*100

/** Severity measure **/
qui reg firstpsyc 1.above c.zabove c.zbelow 1.Dyear (1.above c.zabove c.zbelow)#1.Dyear  [pw=kwt] if  agec!=18 & highrisk==4
global pc_RDDID_highrisk_4=(_b[1.above#1.Dyear]/${overall})*100

qui reg anymental  1.above c.zabove c.zbelow 1.Dyear (1.above c.zabove c.zbelow)#1.Dyear [pw=kwt] if  agec!=18 & highrisk==4
global pt_RDDID_highrisk_4=((_b[_cons]+_b[1.above]+_b[1.Dyear])/${overall_AT})*100


/* Regression Discontinuity */
clear
cd "I:\workdata\706727\Build\Data"
use psychighfreq_AllAges_finalfull

keep if inrange(year,2006,2012)

/* Generate variables */

gen cohab=(!mi(efalle))
gen byte famany_fam=((famany_child>0 & !mi(famany_child)) | (famany_par>0 & !mi(famany_par)) | (famany_sib>0 & !mi(famany_sib)))
gen anymental=(Dpsyc==1 | Dpsychiatry==1 | Dpsykhosp==1 | Dtherapy==1 | Dpsyk_test==1)

/* Merge with own income, education and labor market transfers one year prior */
rename year year2
gen year=year-1
merge m:1 pnr year using incdata, keepusing(pnr year inc pstill isced)
drop if _merge==2
drop _merge

merge m:1 pnr year using J:\Workdata\706727\Temp\Laborcollapse, keepusing(pnr year Dallnonparleave)
drop if _merge==2
drop year _merge
rename year2 year

gen byte EmpEduc=(pstill=="11" | pstill=="12" | pstill=="13" | pstill=="14" | pstill=="19" | pstill=="20" | pstill=="31" | pstill=="32" | pstill=="33" | pstill=="34" | pstill=="35" | pstill=="36"  | pstill=="37" | pstill=="91" )
gen byte highschool=(isced>=3 & !mi(isced))

* Use values from one year prior *
sort pnr year month
by pnr: gen lagfamany=famany[_n-12] if pnr==pnr[_n-12] 
replace famany=lagfamany

by pnr: gen lagfamany_fam=famany_fam[_n-12] if pnr==pnr[_n-12] 
replace famany_fam=lagfamany_fam

by pnr: gen lagcohab=cohab[_n-12] if pnr==pnr[_n-12] 
replace cohab=lagcohab

keep if  !mi(agec) & !mi(koen) & !mi(EmpEduc) & !mi(famany) & !mi(famany_fam) & !mi(inc)  & !mi(Dallnonparleave) & !mi(highschool) & !mi(cohab) & !mi(year)

* Merge with severity prediction *
merge m:1 pnr year using "I:/Workdata/706727/Build/Data/MLSevdata_Pred.dta"
keep if _merge==3
drop _merge

cap drop rankage
gen rankage=floor(agec)

sort year rankage pred
by year rankage: gen double predrank=_n
by year rankage: gen double maxpredrank=_N
replace predrank=ceil((((predrank-0.5)/maxpredrank)*100)/(100/4))
rename predrank highrisk

/* Create income ranks */
set seed 112
gen double random=runiform()
replace random=(random-0.5)/_N

cap drop rankage
gen rankage=floor(agec)

cap drop nonmiss
gen nonmiss=(!mi(inc) & !mi(rankage) & !mi(koen))

sort nonmiss year rankage koen inc random
by nonmiss year rankage koen: gen double incrank=_n
by nonmiss year rankage koen: gen double maxincrank=_N
replace incrank=ceil((((incrank-0.5)/maxincrank)*100+random)/(100/2))
replace incrank=. if nonmiss==0

* Same sample as year_collapse *
drop if year==2012 & month>3 
keep if year>2008 | (year==2008 & month>3)

/* Generate variables */
keep if inrange(agec,20,58)

/* Regression */
g byte above=(agec<38) if !mi(agec)
g byte below=(agec>=38) if !mi(agec)

g double dist=agec-38
g double zabove=above*dist
g double zbelow=dist

g double kwt=max(0,7.5-abs(dist))

/** overall **/
qui reg firstpsyc 1.above c.zabove c.zbelow   [pw=kwt] if  agec!=38
global overall=_b[1.above]

qui reg anymental 1.above c.zabove c.zbelow [pw=kwt] if  agec<38
global overall_AT=_b[_cons]

/** Gender **/
qui reg firstpsyc 1.above c.zabove c.zbelow   [pw=kwt] if  agec!=38 & koen==2
global pc_RD_gender_2=(_b[1.above]/${overall})*100

qui qui reg anymental 1.above c.zabove c.zbelow   [pw=kwt] if  agec!=38 & koen==2
global pt_RD_gender_2=(_b[_cons]/${overall_AT})*100


/** Single **/
qui reg firstpsyc 1.above c.zabove c.zbelow   [pw=kwt] if  agec!=38 & cohab==0
global pc_RD_single=(_b[1.above]/${overall})*100

qui qui reg anymental 1.above c.zabove c.zbelow   [pw=kwt] if  agec!=38 & cohab==0
global pt_RD_single=(_b[_cons]/${overall_AT})*100


/** Income Tertile **/
foreach t in 1  {

qui reg firstpsyc 1.above c.zabove c.zbelow   [pw=kwt] if  agec!=38 & incrank==`t'
global pt_RD_inc_`t'=(_b[1.above]/${overall})*100

qui qui reg anymental 1.above c.zabove c.zbelow   [pw=kwt] if  agec!=38 & incrank==`t'
global pc_RD_inc_`t'=(_b[_cons]/${overall_AT})*100
}

/** High school education **/
qui reg firstpsyc 1.above c.zabove c.zbelow   [pw=kwt] if  agec!=38 & highschool==1
global pc_RD_educ=(_b[1.above]/${overall})*100

qui qui reg anymental 1.above c.zabove c.zbelow   [pw=kwt] if  agec!=38 & highschool==1
global pt_RD_educ=(_b[_cons]/${overall_AT})*100

/** Labor market transfers **/
qui reg firstpsyc 1.above c.zabove c.zbelow   [pw=kwt] if  agec!=38 & Dallnonparleave==1
global pc_RD_labortrans=(_b[1.above]/${overall})*100

qui qui reg anymental 1.above c.zabove c.zbelow   [pw=kwt] if  agec!=38 & Dallnonparleave==1
global pt_RD_labortrans=(_b[_cons]/${overall_AT})*100

/** Own mental health problems **/
qui reg firstpsyc 1.above c.zabove c.zbelow   [pw=kwt] if  agec!=38 & famany==1
global pc_RD_famany=(_b[1.above]/${overall})*100

qui qui reg anymental 1.above c.zabove c.zbelow   [pw=kwt] if  agec!=38 & famany==1
global pt_RD_famany=(_b[_cons]/${overall_AT})*100

/** Own mental health problems **/
qui reg firstpsyc 1.above c.zabove c.zbelow   [pw=kwt] if  agec!=38 & famany_fam==1
global pc_RD_famany_fam=(_b[1.above]/${overall})*100

qui qui reg anymental 1.above c.zabove c.zbelow   [pw=kwt] if  agec!=38 & famany_fam==1
global pt_RD_famany_fam=(_b[_cons]/${overall_AT})*100

/** Severity measure **/
qui reg firstpsyc 1.above c.zabove c.zbelow   [pw=kwt] if  agec!=38 & highrisk==4
global pc_RD_highrisk_4=(_b[1.above]/${overall})*100

qui qui reg anymental 1.above c.zabove c.zbelow   [pw=kwt] if  agec!=38 & highrisk==4
global pt_RD_highrisk_4=(_b[_cons]/${overall_AT})*100


/***** Create Table *******/
clear
mat X=J(8,4,.)

local i=1
foreach var in  gender_2 single inc_1 educ labortrans famany famany_fam highrisk_4 {
	
mat X[`i',1]=${pc_RDDID_`var'}
mat X[`i',2]=${pt_RDDID_`var'}	
mat X[`i',3]=${pc_RD_`var'}
mat X[`i',4]=${pt_RD_`var'}

local i=`i'+1
}

mat Y=X
mat X=X/100

frmttable using I:\workdata\706727\Figs\ComplierTable_CompAll.tex, tex fragment replace statmat(X) sub(1) brackets("","" \ [,]) basefont(fs12) statfont(fs11) coljust(l c) ///
ctitle("Method","RDDID", "RD" \ "Age range", "18", "37")   ///
sdec(2) ///
rtitle("Woman" \" " \ "Single" \ " " \ "Below median income" \" " \  "Upper secondary school or higher" \" " \ "Receive labor market transfers" \" " \"Previous mental illness"  \" " \ "Family mental illness" \" "\ "Top 25 % severity prediction" \" ") ///
multicol(1,2,2)

frmttable using I:\workdata\706727\Figs\ComplierTable_CompAll.rtf, replace statmat(X) sub(1) brackets("","" \ [,] ) basefont(fs12) statfont(fs11)  coljust(l c) ///
ctitle("Method","RDDID", "RD" \ "Age range", "18", "37")   ///
sdec(2) ///
rtitle("Woman" \" " \ "Single" \ " " \ "Below median income" \" " \  "Upper secondary school or higher" \" " \ "Receive labor market transfers" \" " \"Previous mental illness"  \" " \ "Family mental illness" \" " \ "Top 25 % severity prediction" \" ") ///
multicol(1,2,2)
 

/*********************/
/********* 6 *********/
/*********************/

/* RD-DID estimates */
clear
cd "I:\workdata\706727\Build\Data"
use psychighfreq_AllAges_finalfull

keep if inrange(year,2000,2012)

destring(kom),replace
merge m:1 kom year using  "I:\Workdata\706727\Build\Data\ChildhospByKomYear.dta" 
keep if _merge==3
drop _merge

* Same sample as year_collapse *
drop if year==2012 & month>3
keep if inrange(agec,15,23) 

/* Generate variables */
gen AnxDep=(Dbamb_AnxDep==1 | Derp_AnxDep==1 | DPHosp_AnxDep==1)
gen byte firstpsykhospAnxDep=firstpsykhosp*AnxDep==1
gen mort=(!mi(c_liste49))

/* Estimates are not directly comparable to DID results */
foreach var in firstpsyc firstpsykhospAnxDep firstpsychiatry firsttherapy firstpsyk_test firstadep firstanxsleep firstsui mort {
replace `var'=`var'*100*12
}

cap drop Dyear
gen Dyear=0 if inrange(year,2000,2007) | (year==2008 & month<4)
replace Dyear=1 if inrange(year,2009,2012) | (year==2008 & month>=4)

gen cutoff2006=cutoff if year==2006
bysort kom: egen max=max(cutoff2006)
replace cutoff2006=max
drop max

/* Regression */
gen belowcutoff=(agec<cutoff)

g byte above=(agec>=18) if !mi(agec)
g byte below=(agec<18) if !mi(agec)

g double dist=agec-18
g double zabove=above*dist
g double zbelow=dist

g double kwt=max(0,2.5-abs(dist))

foreach var in firstpsyc firstpsykhospAnxDep firstpsychiatry firsttherapy firstpsyk_test firstadep firstanxsleep firstsui mort {

reg `var' 1.above c.zabove c.zbelow 1.Dyear i.cutoff#1.belowcutoff i.cutoff  1.Dyear#(1.above c.zabove c.zbelow)  [pw=kwt] if agec!=18 , vce(cluster pnr)

global est_RDDID_`var'=_b[1.Dyear#1.above]
global se_RDDID_`var'=_se[1.Dyear#1.above]

* Base level *
qui levelsof cutoff, local(lev)

	qui count if e(sample)==1
	local tot=r(N)
	
	global b_RDDID_`var'=_b[_cons]+_b[1.above]+_b[1.Dyear]
	foreach l of local lev {
	qui count if e(sample)==1 & cutoff==`l'
	local share=r(N)

		if (`l'<18) {
		global b_RDDID_`var'=${b_RDDID_`var'}+(`share'/`tot')*_b[`l'.cutoff]
		}
		if (`l'>18) {
		global b_RDDID_`var'=${b_RDDID_`var'}+(`share'/`tot')*(_b[`l'.cutoff]+_b[`l'.cutoff#1.belowcutoff])
		}
	}
}


/* Regression Discontinuity */
clear
cd "I:\workdata\706727\Build\Data"
use psychighfreq_AllAges_finalfull

keep if inrange(year,2000,2012)

* Same sample as year_collapse *
drop if year==2012 & month>3
keep if !inrange(agec,0,18) 

/* Generate variables */
gen AnxDep=(Dbamb_AnxDep==1 | Derp_AnxDep==1 | DPHosp_AnxDep==1)
gen byte firstpsykhospAnxDep=firstpsykhosp*AnxDep==1
gen mort=(!mi(c_liste49))

/* Estimates are not directly comparable to DID results */
foreach var in firstpsyc firstpsykhospAnxDep firstpsychiatry firsttherapy firstpsyk_test firstadep firstanxsleep firstsui mort {
replace `var'=`var'*100*12
}

cap drop Dyear
gen Dyear=0 if inrange(year,2000,2007) | (year==2008 & month<4)
replace Dyear=1 if inrange(year,2009,2012) | (year==2008 & month>=4)

keep if inrange(agec,20,58) & Dyear==1

/* Regression */
g byte above=(agec<38) if !mi(agec)
g byte below=(agec>=38) if !mi(agec)

g double dist=agec-38
g double zabove=above*dist
g double zbelow=dist

g double kwt=max(0,7.5-abs(dist))

foreach var in firstpsyc firstpsykhospAnxDep firstpsychiatry firsttherapy firstpsyk_test firstadep firstanxsleep firstsui mort {

qui reg `var' 1.above c.zabove c.zbelow   [pw=kwt] if Dyear==1 & agec!=38, vce(cluster pnr)

global est_RD_`var'=_b[1.above]
global b_RD_`var'=_b[_cons]
global se_RD_`var'=_se[1.above]
}

/* Create table */
clear

mat X=J(9,6,.)
mat stars=J(9,6,0)

local j=1
foreach met in RDDID RD  {
	
	local i=1
foreach var in firstpsyc firstpsykhospAnxDep firstpsychiatry firsttherapy firstpsyk_test firstadep firstanxsleep firstsui mort {
	
	di "`var'"
	di ${est_`met'_`var'}
	mat X[`i',`j']=${est_`met'_`var'}
	mat X[`i',`j'+1]=${b_`met'_`var'}
	mat X[`i',`j'+2]=${se_`met'_`var'}
	
	mat stars[`i',`j']=(abs(${est_`met'_`var'}/${se_`met'_`var'})>1.645)+(abs(${est_`met'_`var'}/${se_`met'_`var'})>1.96) + (abs(${est_`met'_`var'}/${se_`met'_`var'})>2.58)
	
	local i=`i'+1
	}
local j=`j'+3
}

frmttable using I:\workdata\706727\Figs\Final_FiguresAndTables_MentalHealthTable.tex, tex fragment replace statmat(X) sub(2) annotate(stars) asymbol(*,**,***) brackets("","" \ [,] \ (,)) basefont(fs12) statfont(fs11) coljust(l c) ///
ctitle("Method","RDDID", "RD" \ "Age range", "18",  "37")   ///
sdec(3) ///
rtitle("Psychologist" \" " \ " " \ "Psychiatric hospital" \ " " \ " " \"Psychiatrist"  \" " \ " " \"PCP MHC" \" " \ " " \"PCP Psych. test" \" " \ " " \"Antidepressants" \" " \ " " \"Anxiety and tranquilizers" \" " \ " " \"Suicide attempts" \" " \ " " \"All-cause mortality" \" " \ " ") ///
multicol(1,2,2)

frmttable using I:\workdata\706727\Figs\Final_FiguresAndTables_MentalHealthTable.rtf, replace statmat(X) sub(2) annotate(stars) asymbol(*,**,***) brackets("","" \ [,] \ (,)) basefont(fs12) statfont(fs11)  coljust(l c) ///
ctitle("Method","RDDID", "RD" \ "Age range", "18",  "37")   ///
sdec(3) ///
rtitle("Psychologist" \" " \ " " \ "Psychiatric hospital" \ " " \ " " \"Psychiatrist"  \" " \ " " \"PCP MHC" \" " \ " " \"PCP Psych. test" \" " \ " " \"Antidepressants" \" " \ " " \"Anxiety and tranquilizers" \" " \ " " \"Suicide attempts" \" " \ " " \"All-cause mortality" \" " \ " ") ///
multicol(1,2,2)


/*********************/
/********* 7 *********/
/*********************/

*** Bootstrap standard errors for RD-in-time model ***
clear 
cd "I:\workdata\706727\Build\Data"
use PsykTestLongRunRD, clear

keep if inrange(treatagec,18,50)
keep if inrange(distyrm,-5,7)

keep pnr dupli distyrm treatagec tyrm yrm agec docvisit specdoc medphyc inpat outpat er psyc Dpsyc firstpsyc 

foreach var in psyc docvisit specdoc medphyc inpat outpat er {
replace `var'=0 if `var'<0 | mi(`var')
}

/* Winzorize */
foreach var in psyc docvisit specdoc medphyc inpat outpat er {
qui sum `var' if `var'>0 & !mi(`var'), d
replace `var'=r(p99) if `var'>r(p99) & !mi(`var')
}

egen pnrg=group(pnr dupli)

replace yrm=floor(yrm*4)/4
replace tyrm=floor(tyrm*4)/4
cap drop distyrm*
gen distyrm=yrm-tyrm

gcollapse (mean) psyc firstpsyc docvisit specdoc medphyc inpat outpat er, by(pnr pnrg distyrm treatagec tyrm yrm) fast

g byte treat=(inrange(treatagec,18,37.99))

egen distg=group(distyrm)
egen yrmg=group(yrm)
g byte postreform=(tyrm>=2008.25) & !mi(tyrm)
destring(pnr), gen(numpnr)

/** Create postshock dummies **/

sum distg 
global max=r(max)
global maxm1=r(max)-1

sum distg if distyrm==0
global base=r(max)-1

cap drop distgpoolpre
gen distgpoolpre=distg
replace distgpoolpre=${max} if distg<=${base}
replace distgpoolpre=${base} if distg==${max}

cap drop bgdistg
gen bgdistg=distg
replace bgdistg=${max} if distg==${base}
replace bgdistg=${base} if distg==${max}
replace bgdistg=${max}-1 if distg==1
replace bgdistg=1 if distg==${max}-1

* Find cumulative treatment on individual-level *
gen psyc1=(firstpsyc>0 & !mi(firstpsyc)) if !mi(firstpsyc)
cap drop maxfirstpsyc maxmaxfirstpsyc
bysort pnrg: egen maxfirstpsyc=max(psyc1) if distg>20 & distg<28
by pnrg: egen maxmaxfirstpsyc=max(maxfirstpsyc)

compress

cap program drop aregben Estimate Replay
do "I:\Workdata\706727\STATA Programs\aregben.do"

g byte postshock=(distyrm>=0) & !mi(distyrm)

cap drop d above zabove zbelow 
g double d=tyrm-2008.25
g byte above=(tyrm>=2008.25) if !mi(tyrm)
g double zabove=above*d
g double zbelow=d

cap drop kwt
g double kwt=max(0,7.5-abs(d))

egen tyrmg=group(tyrm)

sum tyrmg if tyrm==2008.25
global tyrm0=r(mean)
sum tyrmg
global tyrmmax=r(max)

keep if inrange(treatagec,18,27.99)
keep if !inrange(tyrm,2008,2008.25)

keep pnr pnrg yrmg kwt bgdistg distgpoolpre zabove zbelow above kwt firstpsyc psyc docvisit specdoc medphyc inpat outpat er
compress

cap postclose john 
postfile john it  str30 var dist koef  using "I:\workdata\706727\Figs\PsykTestLongRun_PooledTable_FullPerDist_RDtime_Bootstrap", replace

forvalues it=1/100 {
preserve
set seed `it'
bsample , cluster(pnr)
	
	foreach var in firstpsyc psyc docvisit specdoc medphyc inpat outpat er   {
		
	qui aregben `var' i.yrmg ib${base}.bgdistg ib${base}.bgdistg#c.zbelow ib${base}.bgdistg#c.zabove ib${base}.distgpoolpre#1.above [pw=kwt], absorb(pnrg) 
	
		forvalues y=20/48 {
		post john (`it')  ("`var'") (`y') (_b[`y'.distgpoolpre#1.above])
		}
	}	
	
di `it'

restore
}
postclose john


/**** Data Preb estimates ****/
clear 
cd "I:\workdata\706727\Build\Data"
use PsykTestLongRunRD, clear

keep if inrange(treatagec,18,50)
keep if inrange(distyrm,-5,7)

keep pnr dupli distyrm treatagec tyrm yrm agec docvisit specdoc medphyc inpat outpat er psyc Dpsyc firstpsyc 

foreach var in psyc docvisit specdoc medphyc inpat outpat er {
replace `var'=0 if `var'<0 | mi(`var')
}

/* Winzorize */
foreach var in psyc docvisit specdoc medphyc inpat outpat er {
qui sum `var' if `var'>0 & !mi(`var'), d
replace `var'=r(p99) if `var'>r(p99) & !mi(`var')
}

egen pnrg=group(pnr dupli)

replace yrm=floor(yrm*4)/4
replace tyrm=floor(tyrm*4)/4
cap drop distyrm*
gen distyrm=yrm-tyrm

gcollapse (mean) psyc firstpsyc docvisit specdoc medphyc inpat outpat er, by(pnr pnrg distyrm treatagec tyrm yrm) fast

g byte treat=(inrange(treatagec,18,37.99))

egen distg=group(distyrm)
egen yrmg=group(yrm)
g byte postreform=(tyrm>=2008.25) & !mi(tyrm)
destring(pnr), gen(numpnr)

/** Create postshock dummies **/

sum distg 
global max=r(max)
global maxm1=r(max)-1

sum distg if distyrm==0
global base=r(max)-1

cap drop distgpoolpre
gen distgpoolpre=distg
replace distgpoolpre=${max} if distg<=${base}
replace distgpoolpre=${base} if distg==${max}

cap drop bgdistg
gen bgdistg=distg
replace bgdistg=${max} if distg==${base}
replace bgdistg=${base} if distg==${max}
replace bgdistg=${max}-1 if distg==1
replace bgdistg=1 if distg==${max}-1

compress

cap program drop aregben Estimate Replay
do "I:\Workdata\706727\STATA Programs\aregben.do"

/**** RD in time ****/

* Standard errors *
preserve
clear 
use "I:\Workdata\706727\Figs\PsykTestLongRun_PooledTable_FullPerDist_RDtime_Bootstrap"

bysort it var: egen sumkoef=sum(koef)

keep it var sumkoef
duplicates drop

gen temp=sumkoef if var=="firstpsyc"
bysort it: egen firststage=max(temp)
	
gen secondstage=sumkoef/firststage

sum sumkoef if  var=="firstpsyc"
global sefirstpsyc=r(sd)

* Second stage *
foreach var in  psyc docvisit specdoc medphyc inpat outpat er   {

sum secondstage if var=="`var'"
global se`var'=r(sd)
}
restore

g byte postshock=(distyrm>=0) & !mi(distyrm)

cap drop d above zabove zbelow 
g double d=tyrm-2008.25
g byte above=(tyrm>=2008.25) if !mi(tyrm)
g double zabove=above*d
g double zbelow=d

cap drop kwt
g double kwt=max(0,7.5-abs(d))

egen tyrmg=group(tyrm)

sum tyrmg if tyrm==2008.25
global tyrm0=r(mean)
sum tyrmg
global tyrmmax=r(max)

mat X=J(8,6,.)
mat stars=J(8,6,0)

preserve
keep if inrange(treatagec,18,27.99)
keep if !inrange(tyrm,2008,2008.25)

* First stage *
local j=1
local i=1

* Baseline 
reg firstpsyc c.zbelow c.zabove 1.above ib${base}.distg ib${base}.distg#c.zbelow ib${base}.distg#c.zabove ib${base}.distg#1.above [pw=kwt]

global baseprob_RDtime=0
forvalues y=20/48 {
global baseprob_RDtime=${baseprob_RDtime} +(_b[_cons]+_b[`y'.distg])*3	
}

* Estimate 
aregben firstpsyc i.yrmg ib${base}.bgdistg ib${base}.bgdistg#c.zbelow ib${base}.bgdistg#c.zabove ib${base}.distgpoolpre#1.above [pw=kwt], absorb(pnrg) vce(cluster pnr)
mat V=e(V)
mat b=e(b)
mat se=J(1,48-20+1,1)*V["20.distgpoolpre#1.above".."48.distgpoolpre#1.above","20.distgpoolpre#1.above".."48.distgpoolpre#1.above"]*J(48-20+1,1,1)
global se=sqrt(se[1,1])

global firststage=0
forvalues y=20/48 {
global firststage=${firststage}+_b[`y'.distgpoolpre#1.above]
}

mat X[`i',`j']=${firststage}*3*100
mat X[`i',`j'+1]=${baseprob_RDtime}*100
mat X[`i',`j'+2]=${se}*3*100
mat stars[`i',`j']=(abs(${firststage}/${se})>1.645)+(abs(${firststage}/${se})>1.96) + (abs(${firststage}/${se})>2.58)

* Baseline 
foreach var in  docvisit specdoc medphyc inpat outpat er  {
		
reg `var' c.zbelow c.zabove 1.above ib${base}.distg ib${base}.distg#c.zbelow ib${base}.distg#c.zabove ib${base}.distg#1.above [pw=kwt]

global baseprob_RDtime_`var'=0
forvalues y=20/48 {
global baseprob_RDtime_`var'=${baseprob_RDtime_`var'} +(_b[_cons]+_b[`y'.distg])*3	
}
}

* Psychologist services *
local j=1
local i=2

aregben psyc i.yrmg ib${base}.bgdistg ib${base}.bgdistg#c.zbelow ib${base}.bgdistg#c.zabove ib${base}.distgpoolpre#1.above [pw=kwt], absorb(pnrg)

global secondstage=0
forvalues y=20/48 {
global secondstage=${secondstage}+_b[`y'.distgpoolpre#1.above]
}

mat X[`i',`j']=${secondstage}/${firststage}
mat X[`i',`j'+2]=${sepsyc}
mat stars[`i',`j']=(abs((${secondstage}/${firststage})/${sepsyc})>1.645)+(abs((${secondstage}/${firststage})/${sepsyc})>1.96) + (abs((${secondstage}/${firststage})/${sepsyc})>2.58)

* Second stage *
local i=3
foreach var in docvisit specdoc medphyc inpat outpat er {

aregben `var' i.yrmg ib${base}.bgdistg ib${base}.bgdistg#c.zbelow ib${base}.bgdistg#c.zabove ib${base}.distgpoolpre#1.above [pw=kwt], absorb(pnrg) 

global secondstage=0
forvalues y=20/48 {
global secondstage=${secondstage}+_b[`y'.distgpoolpre#1.above]
}

mat X[`i',`j']=${secondstage}/${firststage}
mat X[`i',`j'+2]=${se`var'}
mat stars[`i',`j']=(abs((${secondstage}/${firststage})/${se`var'})>1.645)+(abs((${secondstage}/${firststage})/${se`var'})>1.96) + (abs((${secondstage}/${firststage})/${se`var'})>2.58)

mat X[`i',`j'+1]=${baseprob_RDtime_`var'}

local i=`i'+1
}

restore

/** RD 37 **/
keep if inrange(tyrm,2008.25,2012.25)
keep if inrange(treatagec,30,47)

cap drop d above zabove zbelow 
g double d=treatagec-38
g byte above=(treatagec<38) if !mi(treatagec)
g double zabove=above*d
g double zbelow=d

cap drop kwt
g double kwt=max(0,7.5-abs(d))

* First stage *
local j=4
local i=1

* Baseline 
reg firstpsyc c.zbelow c.zabove 1.above ib${base}.distg ib${base}.distg#c.zbelow ib${base}.distg#c.zabove ib${base}.distg#1.above [pw=kwt]

global baseprob_RD=0
forvalues y=20/48 {
global baseprob_RD=${baseprob_RD} +(_b[_cons]+_b[`y'.distg])*3	
}

* Estimate
qui aregben firstpsyc i.yrmg ib${base}.distg ib${base}.distg#c.zbelow ib${base}.distg#c.zabove 1.postshock#1.above  [pw=kwt]  if treatagec!=38, absorb(pnrg) vce(cluster pnr) 

mat X[`i',`j']=_b[1.postshock#1.above]*3*(48-20+1)*100
mat X[`i',`j'+1]=${baseprob_RD}*100
mat X[`i',`j'+2]=_se[1.postshock#1.above]*3*(48-20+1)*100
mat stars[`i',`j']=(abs(_b[1.postshock#1.above]/_se[1.postshock#1.above])>1.645)+(abs(_b[1.postshock#1.above]/_se[1.postshock#1.above])>1.96) + (abs(_b[1.postshock#1.above]/_se[1.postshock#1.above])>2.58)

* Number of psychologist services *
local j=4
local i=2
qui ivreghdfe psyc  ib${base}o1.distg#c.zabove ib${base}o1.distg#c.zbelow (firstpsyc=1.postshock#1.above)  [pw=kwt] if treatage!=38 ,cluster(numpnr)  absorb(pnrg yrmg distg)
eststo reg1

mat X[`i',`j']=_b[firstpsyc]
mat X[`i',`j'+2]=_se[firstpsyc]
mat stars[`i',`j']=(abs(_b[firstpsyc]/_se[firstpsyc])>1.645)+(abs(_b[firstpsyc]/_se[firstpsyc])>1.96) + (abs(_b[firstpsyc]/_se[firstpsyc])>2.58)

* Baseline 
foreach var in  docvisit specdoc medphyc inpat outpat er  {
		
reg `var' c.zbelow c.zabove 1.above ib${base}.distg ib${base}.distg#c.zbelow ib${base}.distg#c.zabove ib${base}.distg#1.above [pw=kwt]

global baseprob_RD_`var'=0
forvalues y=20/48 {
global baseprob_RD_`var'=${baseprob_RD_`var'} +(_b[_cons]+_b[`y'.distg])*3	
}
}

* Second stage *
local i=3
foreach var in docvisit specdoc medphyc inpat outpat er {

qui ivreghdfe `var'  ib${base}o1.distg#c.zabove ib${base}o1.distg#c.zbelow (firstpsyc=1.postshock#1.above)  [pw=kwt] if treatagec!=38 ,cluster(numpnr)  absorb(pnrg yrmg distg)

mat X[`i',`j']=_b[firstpsyc]
mat X[`i',`j'+2]=_se[firstpsyc]
mat stars[`i',`j']=(abs(_b[firstpsyc]/_se[firstpsyc])>1.645)+(abs(_b[firstpsyc]/_se[firstpsyc])>1.96) + (abs(_b[firstpsyc]/_se[firstpsyc])>2.58)

mat X[`i',`j'+1]=${baseprob_RD_`var'}

local i=`i'+1
}

frmttable using "I:\workdata\706727\Figs\LabormarketTable_BothMeth_Final_FullPer.tex", tex fragment replace statmat(X) sub(2) annotate(stars) asymbol(*,**,***) brackets("","" \ [,] \ (,))  coljust(l c)  ///
ctitle("Method","RD in time", "RD" \ "Age range", "18-27", "37")   ///
sdec(2) rtitle("Psychologist (share)" \ "" \ "" \"Psychologist (services)" \ "" \ ""  \ "PCP visit" \ "" \ ""   \  "Specialist visit" \ "" \ "" \  "Somatic drugs" \ "" \ ""  \  "Inpatient visit" \ "" \ "" \  "Outpatient visit" \ "" \ "" \ "ER visit" \ "" \ "" )

frmttable using "I:\workdata\706727\Figs\LabormarketTable_BothMeth_Final_FullPer.rtf",   replace statmat(X) sub(2) annotate(stars) asymbol(*,**,***) brackets("","" \ [,] \ (,))  coljust(l c) ///
ctitle("Method","RD in time", "RD" \ "Age range", "18-27", "37")   ///
sdec(2) rtitle("Psychologist (share)" \ "" \ "" \"Psychologist (services)" \ "" \ ""  \ "PCP visit" \ "" \ ""   \  "Specialist visit" \ "" \ "" \  "Somatic drugs" \ "" \ ""  \  "Inpatient visit" \ "" \ "" \  "Outpatient visit" \ "" \ "" \ "ER visit" \ "" \ "" )


/*********************/
/********* 8 *********/
/*********************/

/*** Bootstrap standard errors for RD-in-time model ***/
clear all
cd "I:\workdata\706727\Build\Data"
use PsykTestLongRunRD, clear

keep if inrange(treatagec,21,27.99)
keep if inrange(distyrm,-3.5,7)

** drop individuals just before and at cutoff ** 
drop if inrange(tyrm,2008,2008.25)

keep pnr dupli distyrm treatagec tyrm yrm agec firstpsyc Dpsyc Ddisab Dsick Dsueduc Dallnonparleave emp

drop if agec<=18

egen pnrg=group(pnr dupli)

replace yrm=floor(yrm*4)/4
gen tyrm2=floor(tyrm*4)/4
cap drop distyrm*
gen distyrm=yrm-tyrm2

keep if inrange(distyrm,-2.75,7)
gcollapse (mean) firstpsyc Dpsyc Ddisab Dsick Dsueduc Dallnonparleave emp, by(pnr pnrg distyrm treatagec tyrm yrm) fast

egen distg=group(distyrm)
egen yrmg=group(yrm)

/** Create postshock dummies **/
sum distg 
global max=r(max)
global maxm1=r(max)-1

/* set base group to three periods before treatment */
sum distg if distyrm==0
global bdist=1
global base=r(max)-${bdist}

cap drop distgpoolpre
gen distgpoolpre=distg
replace distgpoolpre=${max} if distg<=${base}
replace distgpoolpre=${base} if distg==${max}

cap drop bgdistg
gen bgdistg=distg
replace bgdistg=${max} if distg==${base}
replace bgdistg=${base} if distg==${max}
replace bgdistg=${max}-1 if distg==1
replace bgdistg=1 if distg==${max}-1

cap program drop aregben Estimate Replay
do "I:\Workdata\706727\STATA Programs\aregben.do"

compress

/**** RD in time ****/
g byte postshock=(distyrm>=0) & !mi(distyrm)

cap drop d above zabove zbelow 
g double d=tyrm-2008.25
g byte above=(tyrm>=2008.25) if !mi(tyrm)
g double zabove=above*d
g double zbelow=d

cap drop kwt
g double kwt=max(0,7.5-abs(d))

gen tyrm2=floor(tyrm*4)/4

keep pnr pnrg yrmg tyrm2 kwt bgdistg distgpoolpre zabove zbelow above kwt firstpsyc Ddisab Dsick Dsueduc Dallnonparleave emp
compress

cap postclose john 
postfile john it  str30 var dist koef  using "I:\workdata\706727\Figs\PsykTestLongRun_MonthlyLabor_Final_FullPer_RDtime_Bootstrap", replace

forvalues it=1/100 {
preserve
set seed `it'
bsample , cluster(pnr) 
	
	foreach var in firstpsyc Ddisab Dsick Dsueduc Dallnonparleave emp   {
		
	qui aregben `var' i.yrmg ib${base}.bgdistg ib${base}.bgdistg#c.zbelow ib${base}.bgdistg#c.zabove ib${base}.distgpoolpre#1.above [pw=kwt], absorb(pnrg) 
	
		forvalues y=11/39 {
		post john (`it')  ("`var'") (`y') (_b[`y'.distgpoolpre#1.above])
		}
	}	
	
di `it'

restore
}
postclose john

/**** RD-in-Time Estimates ****/
clear all
cd "I:\workdata\706727\Build\Data"
use PsykTestLongRunRD, clear

keep if inrange(treatagec,21,27.99)
keep if inrange(distyrm,-3.5,7)

** drop individuals just before and at cutoff ** 
drop if inrange(tyrm,2008,2008.25)

keep pnr dupli distyrm treatagec tyrm yrm agec firstpsyc Dpsyc Ddisab Dsick Dsueduc Dallnonparleave emp

drop if agec<=18

egen pnrg=group(pnr dupli)

replace yrm=floor(yrm*4)/4
gen tyrm2=floor(tyrm*4)/4
cap drop distyrm*
gen distyrm=yrm-tyrm2

keep if inrange(distyrm,-2.75,7)
gcollapse (mean) firstpsyc Dpsyc Ddisab Dsick Dsueduc Dallnonparleave emp, by(pnr pnrg distyrm treatagec tyrm yrm) fast

g byte postshock=(distyrm>=0) & !mi(distyrm)

cap drop d above zabove zbelow 
g double d=tyrm-2008.25
g byte above=(tyrm>=2008.25) if !mi(tyrm)
g double zabove=above*d
g double zbelow=d

cap drop kwt
g double kwt=max(0,7.5-abs(d))

egen distg=group(distyrm)
egen yrmg=group(yrm)
egen tyrmg=group(tyrm)

sum distg 
global max=r(max)
global maxm1=r(max)-1

/* set base group to three periods before treatment */
sum distg if distyrm==0
global bdist=1
global base=r(max)-${bdist}

cap drop distgpoolpre
gen distgpoolpre=distg
replace distgpoolpre=${max} if distg<=${base}
replace distgpoolpre=${base} if distg==${max}

cap drop bgdistg
gen bgdistg=distg
replace bgdistg=${max} if distg==${base}
replace bgdistg=${base} if distg==${max}
replace bgdistg=${max}-1 if distg==1
replace bgdistg=1 if distg==${max}-1

* Standard errors *
preserve
clear 
use "I:\Workdata\706727\Figs\PsykTestLongRun_MonthlyLabor_Final_FullPer_RDtime_Bootstrap"

bysort it var: egen sumkoef=sum(koef)

keep it var sumkoef
duplicates drop

gen temp=sumkoef if var=="firstpsyc"
bysort it: egen firststage=max(temp)
	
gen secondstage=sumkoef/firststage

sum sumkoef if  var=="firstpsyc"
global sefirstpsyc=r(sd)

* Second stage *
foreach var in  Ddisab Dsick Dsueduc Dallnonparleave emp  {

sum secondstage if var=="`var'"
global se`var'=r(sd)
}
restore

cap program drop aregben Estimate Replay
do "I:\Workdata\706727\STATA Programs\aregben.do"


* First stage *
local i=1
local j=1

aregben firstpsyc i.yrmg ib${base}.bgdistg ib${base}.bgdistg#c.zbelow ib${base}.bgdistg#c.zabove ib${max}.distgpoolpre#1.above [pw=kwt], absorb(pnrg) vce(cluster pnr)
mat V=e(V)
mat b=e(b)
mat se=J(1,39-11+1,1)*V["11.distgpoolpre#1.above".."39.distgpoolpre#1.above","11.distgpoolpre#1.above".."39.distgpoolpre#1.above"]*J(39-11+1,1,1)
mat est=b[1,"11.distgpoolpre#1.above".."39.distgpoolpre#1.above"]*J(39-11+1,1,1)

global firststage=est[1,1]
global se_RDtime_firstpsyc=sqrt(se[1,1])*3*100
global est_RDtime_firstpsyc=b[1,1]*3*100

* Baseline *
foreach var in Dsick Ddisab Dsueduc Dallnonparleave emp {
		
reg `var' c.zbelow c.zabove 1.above ib${base}.distg ib${base}.distg#c.zbelow ib${base}.distg#c.zabove ib${base}.distg#1.above [pw=kwt]

global b_RDtime_`var'=0
forvalues y=11/39 {
global b_RDtime_`var'=${b_RDtime_`var'} +(_b[_cons]+_b[`y'.distg])*3	
}
}

* Second stage *
destring(pnr), gen(numpnr)

foreach var in  Dsick Ddisab Dsueduc Dallnonparleave emp {
	
qui aregben `var'  i.yrmg ib${base}.bgdistg ib${base}.bgdistg#c.zbelow ib${base}.bgdistg#c.zabove ib${max}.distgpoolpre#1.above [pw=kwt], absorb(pnrg) vce(cluster pnr)
mat b=e(b)
mat est=b[1,"11.distgpoolpre#1.above".."39.distgpoolpre#1.above"]*J(39-11+1,1,1)

global est_RDtime_`var'=est[1,1]/${firststage}
global se_RDtime_`var'=${se`var'}

}

/** RD 37 **/
clear 
cd "I:\workdata\706727\Build\Data"
use PsykTestLongRunRD, clear

keep if inrange(tyrm,2008.25,2012.25)
keep if inrange(treatagec,30,47)
keep if inrange(distyrm,-5,7)

keep pnr dupli distyrm treatagec tyrm yrm agec firstpsyc Dpsyc Ddisab Dsick Dsueduc Dallnonparleave emp

egen pnrg=group(pnr dupli)

replace yrm=floor(yrm*4)/4
gen tyrm2=floor(tyrm*4)/4
cap drop distyrm*
gen distyrm=yrm-tyrm2

gcollapse (mean) firstpsyc Dpsyc Ddisab Dsick Dsueduc Dallnonparleave emp, by(pnr pnrg distyrm treatagec tyrm yrm) fast

g byte postshock=(distyrm>=0) & !mi(distyrm)
egen distg=group(distyrm)
egen yrmg=group(yrm)

cap drop d above zabove zbelow 
g double d=treatagec-38
g byte above=(treatagec<38) if !mi(treatagec)
g double zabove=above*d
g double zbelow=d

cap drop kwt
g double kwt=max(0,7.5-abs(d))

* Baseline *
foreach var in Dsick Ddisab Dallnonparleave emp {
		
reg `var' c.zbelow c.zabove 1.above ib${base}.distg ib${base}.distg#c.zbelow ib${base}.distg#c.zabove ib${base}.distg#1.above [pw=kwt]

global b_RD_`var'=0
forvalues y=20/48 {
global b_RD_`var'=${b_RD_`var'} +(_b[_cons]+_b[`y'.distg])*3	
}
}

* Second stage *
destring(pnr), gen(numpnr)

foreach var in  Dsick Ddisab Dallnonparleave emp {

qui ivreghdfe `var'  ib${base}o1.distg#c.zabove ib${base}o1.distg#c.zbelow (firstpsyc=1.postshock#1.above)  [pw=kwt] if treatage!=38 ,cluster(numpnr)  absorb(pnrg yrmg distg)

global est_RD_`var'=_b[firstpsyc]
global se_RD_`var'=_se[firstpsyc]

}

/**** Create Table ****/
clear

mat X=J(5,6,.)
mat stars=J(5,6,0)

local j=1
foreach met in RDtime RD  {
	
	local i=1
foreach var in  Dsick Ddisab Dsueduc Dallnonparleave emp {
	
	cap mat X[`i',`j']=${est_`met'_`var'}
	cap mat X[`i',`j'+1]=${b_`met'_`var'}
	cap mat X[`i',`j'+2]=${se_`met'_`var'}

	cap mat stars[`i',`j']=(abs(${est_`met'_`var'}/${se_`met'_`var'})>1.645)+(abs(${est_`met'_`var'}/${se_`met'_`var'})>1.96) + (abs(${est_`met'_`var'}/${se_`met'_`var'})>2.58)
	
	local i=`i'+1
	}
local j=`j'+3
}

frmttable using "I:\workdata\706727\Figs\LabormarketTable_MonthlyLaborOutcomes_Final_FullPer_new.tex", tex fragment replace statmat(X) sub(2) annotate(stars) asymbol(*,**,***) brackets("","" \ [,] \ (,))  coljust(l c)  ///
ctitle("Method","RD in time", "RD" \ "Age range", "21-27", "37")   ///
sdec(2) rtitle( "Sick leave"\ "" \ "" \ "Disability pension" \ "" \ "" \ "Under education" \ "" \ "" \ "Any labor market transfer" \ " " \ " " \ "Employment" \ "" \ "")

frmttable using "I:\workdata\706727\Figs\LabormarketTable_MonthlyLaborOutcomes_Final_FullPer_new.rtf",   replace statmat(X) sub(2) annotate(stars) asymbol(*,**,***) brackets("","" \ [,] \ (,))  coljust(l c) ///
ctitle("Method","RD in time", "RD" \ "Age range", "21-27", "37")   ///
sdec(2) rtitle( "Sick leave"\ "" \ "" \ "Disability pension" \ "" \ ""  \ "Under education" \ "" \ "" \ "Any labor market transfer" \ " " \ " " \ "Employment" \ "" \ "")


/*********************/
/********* 9 *********/
/*********************/

/*** Bootstrap standard errors for RD-in-time model ***/
clear 
cd "I:\workdata\706727\Build\Data"
use PsykTestLongRunRD, clear

keep if inrange(treatagec,18,50)
keep if inrange(distyrm,-5,7)

keep pnr dupli distyrm treatagec tyrm year month yrm agec docvisit Ddocvisit p_docvisit inpat Dinpat p_inpat outpat Doutpat p_outpat er Der p_er medphyc Dmedphyc p_medphyc specdoc Dspecdoc p_specdoc psyc p_psyc Dpsyc firstpsyc 

foreach var in psyc docvisit specdoc medphyc inpat outpat er {
replace `var'=0 if mi(`var') | `var'<0
replace p_`var'=0 if mi(p_`var') | p_`var'<0
}

gen year2=year if inrange(month,4,12)
replace year2=year-1 if inrange(month,1,3)

foreach var in p_psyc p_docvisit p_specdoc  {
replace `var'=`var'/100 if year<=2004
}
** 2020 prices in primary care (based on reimbursement for GP consultation) **
foreach var in p_psyc p_docvisit p_specdoc p_medphyc {

replace `var'=`var'*(145.46/97.11) if year2==2000
replace `var'=`var'*(145.46/100.35) if year2==2001
replace `var'=`var'*(145.46/102.59) if year2==2002
replace `var'=`var'*(145.46/105.57) if year2==2003
replace `var'=`var'*(145.46/106.81) if year2==2004
replace `var'=`var'*(145.46/108.05) if year2==2005
replace `var'=`var'*(145.46/113.25) if year2==2006
replace `var'=`var'*(145.46/116.90) if year2==2007
replace `var'=`var'*(145.46/123.02) if year2==2008
replace `var'=`var'*(145.46/126.86) if year2==2009
replace `var'=`var'*(145.46/128.76) if year2==2010
replace `var'=`var'*(145.46/129.52) if year2==2011
replace `var'=`var'*(145.46/131.93) if year2==2012
replace `var'=`var'*(145.46/132.44) if year2==2013
replace `var'=`var'*(145.46/133.84) if year2==2014
replace `var'=`var'*(145.46/134.92) if year2==2015
replace `var'=`var'*(145.46/136.66) if year2==2016
replace `var'=`var'*(145.46/138.94) if year2==2017
replace `var'=`var'*(145.46/140.67) if year2==2018
replace `var'=`var'*(145.46/142.10) if year2==2019
}

* Winzorize *
foreach var in psyc docvisit specdoc medphyc inpat outpat er {
qui sum p_`var' if p_`var'>0 & !mi(p_`var') , d
replace p_`var'=r(p99) if p_`var'>r(p99) & !mi(p_`var')
}

keep if inrange(treatagec,18,27.99)
keep if !inrange(tyrm,2008,2008.25)

gen p_allphyc=p_docvisit + p_specdoc + p_medphyc + p_inpat + p_outpat + p_er

egen pnrg=group(pnr dupli)

replace yrm=floor(yrm*4)/4
replace tyrm=floor(tyrm*4)/4
cap drop distyrm*
gen distyrm=yrm-tyrm

gcollapse (mean)  firstpsyc p_psyc p_inpat p_outpat p_er p_allphyc, by(pnr pnrg distyrm treatagec tyrm yrm) fast

egen distg=group(distyrm)
egen yrmg=group(yrm)

/** Discount **/
global ry=0.035
global rq=(1+${ry})^(1/4) -1
gen quart=floor(distyrm*4)

foreach var in psyc inpat outpat er allphyc {
replace p_`var'=p_`var'/((1+${rq})^(quart)) if quart>0 & !mi(quart)
}

/**** RD in time ****/
g byte postshock=(distyrm>=0) & !mi(distyrm)

cap drop d above zabove zbelow 
g double d=tyrm-2008.25
g byte above=(tyrm>=2008.25) if !mi(tyrm)
g double zabove=above*d
g double zbelow=d

cap drop kwt
g double kwt=max(0,7.5-abs(d))

egen tyrmg=group(tyrm)

sum tyrmg if tyrm==2008.25
global tyrm0=r(mean)
sum tyrmg
global tyrmmax=r(max)


/** Create postshock dummies **/
sum distg 
global max=r(max)
global maxm1=r(max)-1

sum distg if distyrm==0
global base=r(max)-1

cap drop distgpoolpre
gen distgpoolpre=distg
replace distgpoolpre=${max} if distg<=${base}
replace distgpoolpre=${base} if distg==${max}

cap drop bgdistg
gen bgdistg=distg
replace bgdistg=${max} if distg==${base}
replace bgdistg=${base} if distg==${max}
replace bgdistg=${max}-1 if distg==1
replace bgdistg=1 if distg==${max}-1

keep pnr pnrg yrmg kwt bgdistg distgpoolpre zabove zbelow above kwt firstpsyc p_psyc p_inpat p_outpat p_er p_allphyc
compress

cap program drop aregben Estimate Replay
do "I:\Workdata\706727\STATA Programs\aregben.do"

cap postclose john 
postfile john it  str30 var dist koef  using "I:\workdata\706727\Figs\PsykTestLongRun_CostTable_FullPerDist_RDtime_Bootstrap", replace

forvalues it=1/100 {
preserve
set seed `it'
bsample , cluster(pnr)
	
	foreach var in firstpsyc p_psyc p_inpat p_outpat p_er p_allphyc  {
		
	qui aregben `var' i.yrmg ib${base}.bgdistg ib${base}.bgdistg#c.zbelow ib${base}.bgdistg#c.zabove ib${base}.distgpoolpre#1.above [pw=kwt], absorb(pnrg) 
	
		forvalues y=20/48 {
		post john (`it')  ("`var'") (`y') (_b[`y'.distgpoolpre#1.above])
		}
	}	
	
di `it'

restore
}
postclose john


/**** Data Preb estimates ****/
clear 
cd "I:\workdata\706727\Build\Data"
use PsykTestLongRunRD, clear

keep if inrange(treatagec,18,50)
keep if inrange(distyrm,-5,7)

keep pnr dupli distyrm treatagec tyrm year month yrm agec docvisit Ddocvisit p_docvisit inpat Dinpat p_inpat outpat Doutpat p_outpat er Der p_er medphyc Dmedphyc p_medphyc specdoc Dspecdoc p_specdoc psyc p_psyc Dpsyc firstpsyc 

foreach var in psyc docvisit specdoc medphyc inpat outpat er {
replace `var'=0 if mi(`var') | `var'<0
replace p_`var'=0 if mi(p_`var') | p_`var'<0
}

gen year2=year if inrange(month,4,12)
replace year2=year-1 if inrange(month,1,3)

foreach var in p_psyc p_docvisit p_specdoc  {
replace `var'=`var'/100 if year<=2004
}
** 2020 prices in primary care (based on reimbursement for GP consultation) **
foreach var in p_psyc p_docvisit p_specdoc p_medphyc {

replace `var'=`var'*(145.46/97.11) if year2==2000
replace `var'=`var'*(145.46/100.35) if year2==2001
replace `var'=`var'*(145.46/102.59) if year2==2002
replace `var'=`var'*(145.46/105.57) if year2==2003
replace `var'=`var'*(145.46/106.81) if year2==2004
replace `var'=`var'*(145.46/108.05) if year2==2005
replace `var'=`var'*(145.46/113.25) if year2==2006
replace `var'=`var'*(145.46/116.90) if year2==2007
replace `var'=`var'*(145.46/123.02) if year2==2008
replace `var'=`var'*(145.46/126.86) if year2==2009
replace `var'=`var'*(145.46/128.76) if year2==2010
replace `var'=`var'*(145.46/129.52) if year2==2011
replace `var'=`var'*(145.46/131.93) if year2==2012
replace `var'=`var'*(145.46/132.44) if year2==2013
replace `var'=`var'*(145.46/133.84) if year2==2014
replace `var'=`var'*(145.46/134.92) if year2==2015
replace `var'=`var'*(145.46/136.66) if year2==2016
replace `var'=`var'*(145.46/138.94) if year2==2017
replace `var'=`var'*(145.46/140.67) if year2==2018
replace `var'=`var'*(145.46/142.10) if year2==2019
}

* Winzorize *
foreach var in psyc docvisit specdoc medphyc inpat outpat er {
qui sum p_`var' if p_`var'>0 & !mi(p_`var') , d
replace p_`var'=r(p99) if p_`var'>r(p99) & !mi(p_`var')
}

gen p_allphyc=p_docvisit + p_specdoc + p_medphyc + p_inpat + p_outpat + p_er

egen pnrg=group(pnr dupli)

replace yrm=floor(yrm*4)/4
replace tyrm=floor(tyrm*4)/4
cap drop distyrm*
gen distyrm=yrm-tyrm

gcollapse (mean)  firstpsyc p_psyc p_docvisit p_specdoc p_medphyc p_inpat p_outpat p_er p_allphyc, by(pnr pnrg distyrm treatagec tyrm yrm) fast

/**** RD in time ****/

* Standard errors *
preserve
clear 
use "I:\Workdata\706727\Figs\PsykTestLongRun_CostTable_FullPerDist_RDtime_Bootstrap"

bysort it var: egen sumkoef=sum(koef)

keep it var sumkoef
duplicates drop

gen temp=sumkoef if var=="firstpsyc"
bysort it: egen firststage=max(temp)

cap drop temp
gen temp=sumkoef if var=="p_psyc"
bysort it: egen firststage2=max(temp)
	
gen secondstage=sumkoef/firststage
gen offset=sumkoef/firststage2

* Second stage *
foreach var in  psyc inpat outpat er   {

sum secondstage if var=="p_`var'"
global se`var'=r(sd)
}
sum offset if var=="p_allphyc"
global seoffset=r(sd)
restore

g byte postshock=(distyrm>=0) & !mi(distyrm)

cap drop d above zabove zbelow 
g double d=tyrm-2008.25
g byte above=(tyrm>=2008.25) if !mi(tyrm)
g double zabove=above*d
g double zbelow=d

cap drop kwt
g double kwt=max(0,7.5-abs(d))

egen distg=group(distyrm)
egen yrmg=group(yrm)
egen tyrmg=group(tyrm)

sum tyrmg if tyrm==2008.25
global tyrm0=r(mean)
sum tyrmg
global tyrmmax=r(max)

/** Create postshock dummies **/
sum distg 
global max=r(max)
global maxm1=r(max)-1

sum distg if distyrm==0
global base=r(max)-1

cap drop distgpoolpre
gen distgpoolpre=distg
replace distgpoolpre=${max} if distg<=${base}
replace distgpoolpre=${base} if distg==${max}

cap drop bgdistg
gen bgdistg=distg
replace bgdistg=${max} if distg==${base}
replace bgdistg=${base} if distg==${max}
replace bgdistg=${max}-1 if distg==1
replace bgdistg=1 if distg==${max}-1

compress

/** Discount **/
global ry=0.035
global rq=(1+${ry})^(1/4) -1
gen quart=floor(distyrm*4)

foreach var in psyc docvisit specdoc medphyc inpat outpat er allphyc {
replace p_`var'=p_`var'/((1+${rq})^(quart)) if quart>0 & !mi(quart)
}

mat X=J(5,4,.)
mat stars=J(5,4,0)

cap program drop aregben Estimate Replay
do "I:\Workdata\706727\STATA Programs\aregben.do"

preserve
keep if inrange(treatagec,18,27.99)

** drop individuals just before cutoff ** 
drop if inrange(tyrm,2008,2008.25)

** First stage 1 **
qui aregben firstpsyc i.yrmg ib${base}.bgdistg ib${base}.bgdistg#c.zbelow ib${base}.bgdistg#c.zabove ib${base}.distgpoolpre#1.above [pw=kwt], absorb(pnrg) 
mat b=e(b)
mat est=b[1,"20.distgpoolpre#1.above".."48.distgpoolpre#1.above"]*J(48-20+1,1,1)
global firststage=est[1,1]

** First stage 2 **
qui aregben p_psyc i.yrmg ib${base}.bgdistg ib${base}.bgdistg#c.zbelow ib${base}.bgdistg#c.zabove ib${base}.distgpoolpre#1.above [pw=kwt], absorb(pnrg)  
mat b=e(b)
mat est=b[1,"20.distgpoolpre#1.above".."48.distgpoolpre#1.above"]*J(48-20+1,1,1)
global firststage2=est[1,1]

local i=1
local j=1
mat X[`i',`j']=est[1,1]/${firststage}
mat X[`i',`j'+1]=${sepsyc}
local ttest=(est[1,1]/${firststage})/${sepsyc}
mat stars[`i',`j']=(abs(`ttest')>1.645)+(abs(`ttest')>1.96) + (abs(`ttest')>2.58)

* Remaining outcomes *
local i=2
local j=1
foreach var in inpat outpat er  {

qui aregben p_`var' i.yrmg ib${base}.bgdistg ib${base}.bgdistg#c.zbelow ib${base}.bgdistg#c.zabove ib${base}.distgpoolpre#1.above [pw=kwt], absorb(pnrg) 
mat b=e(b)
mat est=b[1,"20.distgpoolpre#1.above".."48.distgpoolpre#1.above"]*J(48-20+1,1,1)

mat X[`i',`j']=est[1,1]/${firststage}
mat X[`i',`j'+1]=${se`var'}
local ttest=(est[1,1]/${firststage})/${se`var'}
mat stars[`i',`j']=(abs(`ttest')>1.645)+(abs(`ttest')>1.96) + (abs(`ttest')>2.58)

local i=`i'+1
}

* Offset *
qui aregben p_allphyc i.yrmg ib${base}.bgdistg ib${base}.bgdistg#c.zbelow ib${base}.bgdistg#c.zabove ib${base}.distgpoolpre#1.above [pw=kwt], absorb(pnrg)  
mat b=e(b)
mat est=b[1,"20.distgpoolpre#1.above".."48.distgpoolpre#1.above"]*J(48-20+1,1,1)

mat X[`i',`j']=est[1,1]/${firststage2}
mat X[`i',`j'+1]=${seoffset}
local ttest=(est[1,1]/${firststage2})/${seoffset}
mat stars[`i',`j']=(abs(`ttest')>1.645)+(abs(`ttest')>1.96) + (abs(`ttest')>2.58)

restore

/** RD 37 **/
keep if inrange(tyrm,2008.25,2012.25)
keep if inrange(treatagec,30,47)

cap drop above zabove zbelow 
cap drop d
g double d=treatagec-38
g byte above=(treatagec<38) if !mi(treatagec)
g double zabove=above*d
g double zbelow=d

cap drop kwt
g double kwt=max(0,7.5-abs(d))

* First stage *
local j=3
local i=1
qui ivreghdfe p_psyc  i.yrmg#c.zabove i.yrmg#c.zbelow  ib${base}.distg#c.zabove ib${base}.distg#c.zbelow (firstpsyc=1.postshock#1.above) [pw=kwt] if treatage!=38, cluster(pnr) absorb(pnrg yrmg distg)

mat X[`i',`j']=_b[firstpsyc]
mat X[`i',`j'+1]=_se[firstpsyc]
mat stars[`i',`j']=(abs(_b[firstpsyc]/_se[firstpsyc])>1.645)+(abs(_b[firstpsyc]/_se[firstpsyc])>1.96) + (abs(_b[firstpsyc]/_se[firstpsyc])>2.58)

* Other variables *
local i=2
foreach var in  inpat outpat er   {

 ivreghdfe p_`var'  i.yrmg#c.zabove i.yrmg#c.zbelow  ib${base}o1.distg#c.zabove ib${base}o1.distg#c.zbelow (firstpsyc=1.postshock#1.above) [pw=kwt] if treatage!=38, cluster(pnr) absorb(pnrg yrmg distg)

mat X[`i',`j']=_b[firstpsyc]
mat X[`i',`j'+1]=_se[firstpsyc]
mat stars[`i',`j']=(abs(_b[firstpsyc]/_se[firstpsyc])>1.645)+(abs(_b[firstpsyc]/_se[firstpsyc])>1.96) + (abs(_b[firstpsyc]/_se[firstpsyc])>2.58)

local i=`i'+1
}

* Offset *
qui ivreghdfe p_allphyc  i.yrmg#c.zabove i.yrmg#c.zbelow  ib${base}o1.distg#c.zabove ib${base}o1.distg#c.zbelow (p_psyc=1.postshock#1.above) [pw=kwt] if treatage!=38, cluster(pnr) absorb(pnrg yrmg distg)

mat X[`i',`j']=_b[p_psyc]
mat X[`i',`j'+1]=_se[p_psyc]
mat stars[`i',`j']=(abs(_b[p_psyc]/_se[p_psyc])>1.645)+(abs(_b[p_psyc]/_se[p_psyc])>1.96) + (abs(_b[p_psyc]/_se[p_psyc])>2.58)


* Create table *
frmttable using "I:\workdata\706727\Figs\LabormarketTable_New_BothMeth_CostTable_New_FullPer.tex", tex fragment replace statmat(X) sub(1) annotate(stars) asymbol(*,**,***)  brackets("","" \ (,))  coljust(l c)  ///
ctitle("Method","RD in time", "RD" \ "Age range", "18-27", "37")   ///
sdec(0\0\0\0\0\0\0\0\2\2) rtitle("Psychologist" \ ""   \ "Inpatient visit" \ "" \  "Outpatient visit" \ "" \ "ER visit" \ "" \ "Physical health offset"\ " " )

frmttable using "I:\workdata\706727\Figs\LabormarketTable_New_BothMeth_CostTable_New_FullPer.rtf",   replace statmat(X) sub(1) annotate(stars) asymbol(*,**,***)  brackets("","" \ (,)) coljust(l c) ///
ctitle("Method","RD in time", "RD" \ "Age range", "18-27", "37")   ///
sdec(0\0\0\0\0\0\0\0\2\2) rtitle("Psychologist" \ ""   \ "Inpatient visit" \ "" \  "Outpatient visit" \ "" \ "ER visit" \ "" \ "Physical health offset"\ " " )








